1. Symmetric Stiffness Tensor and Voigt notation

#fem/stiffness-tensor-symmetry#fem/symmetric-index-mapping#fem/voigt-notation

1. Concepts

(1) Symmetric relations in

Voigt notation is a simple way to represent a symmetric tensor.
For a $2\times 2$ symmetric matrix with $x_{12}=x_{21}$:

$$\begin{array}{}\text{(1.1.1)}& X=\left[\begin{array}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\end{array}\right]\stackrel{\text{voigt index}}{⟶}\left[\begin{array}{c}{x}_{11}\\ x_{22}\\ x_{12}\end{array}\right]\end{array}$$

0 -> (0,0)
1 -> (1,1)
2 -> (0,1) and (1,0)

For $3\times 3$ symmetric strain matrix. There's 6 independent coordinates :

$$\begin{array}{}\text{(1.1.2)}& \sigma =\left[\begin{array}{ccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}\\ {\sigma }_{21}& {\sigma }_{22}& {\sigma }_{23}\\ {\sigma }_{31}& {\sigma }_{32}& {\sigma }_{33}\end{array}\right]=\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]\text{where}{\sigma }_4={\sigma }_{23}{\sigma }_5={\sigma }_{13}{\sigma }_6={\sigma }_{12}\end{array}$$

0 -> (0,0)
1 -> (1,1)
2 -> (2,2)
3 -> (1,2) and (2,1) 
4 -> (0,2) and (2,0) 
5 -> (0,1) and (1,0)

Conventional sequence We often use ${\sigma }_4={\sigma }_{23}{\sigma }_5={\sigma }_{13}{\sigma }_6={\sigma }_{12}$ here.

For general formula for the number of Voigt indices can be calculated by

$$\begin{array}{}\text{(1.1.3)}& \text{Number of Voigt indices}=\left(\begin{array}{c}d+n-1\\ n\end{array}\right)=C_{d+n-1}^n\end{array}$$

Note here we know from

$$\begin{array}{}\text{(1.1.4)}& A(n,r)=\frac{n!}{(n-r)!}C(n,r)=\frac{n!}{(n-r)!r!}\end{array}$$

where tensor rank is $n$ and $d$ is dimension of space.

(2) Symmetry of stiffness tensor

For elastic tensors $E_{ijkl}$ It's a 4-th order tensor, which has 3x3x3x3 = 81 elements. In a general case, $E_{ijkl}$ is generally do not have this full symmetry.But under the assumption of symmetric stress and strain :

$$\begin{array}{}\text{(1.2.1)}& {\sigma }_{ij}={\sigma }_{ji}{\epsilon }_{ij}={\epsilon }_{ji}\end{array}$$

We have 2 minor symmetries :

$$\begin{array}{}\text{(1.2.2)}& C_{ijkl}=C_{jikl}{C}_{ijkl}=C_{ijlk}\end{array}$$

And also under a strain-energy relation :

$$\begin{array}{}\text{(1.2.3)}& W=\frac{1}{2}{C}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}=\frac{1}{2}{C}_{klij}{\epsilon }_{kl}{\epsilon }_{ij}{C}_{ijkl}=\frac{{\partial }^{2}W}{\partial {\epsilon }_{ij}\partial {\epsilon }_{kl}}=C_{klij}\end{array}$$

So a general anisotropic elastic material usually satisfy the minor and major symmetries :

$$\begin{array}{}\text{(1.2.4)}& \boxed{{\displaystyle C_{ijkl}=C_{jikl}=C_{ijlk}=C_{klij}}}\end{array}$$

We can use a symmetric $6\times 6$ matrix with $21$ independent constants to express the stiffness of a general anisotropic elastic material.

So by setting $ij=I,kl=J$, We can also use Voigt Index to map the $9\times 9$ stiffness matrix into a $6\times 6$ one, as shown in Fig. 1.1. We also note, in the composite mechanics, the general Hook's Law ^[1] use a 6x6 matrix.

6x6 stiffness matrix in voigt notation

which is 3x3 matrix (6 independent variable)

The sequence in picture is not same as conventional sequence, but valid for a

$${\sigma }_4={\sigma }_{01}{\sigma }_5={\sigma }_{02}{\sigma }_6={\sigma }_{12}$$

which is the sequence we later applied in the programming.

Why we often use the conventional sequence is because the consistency for 2D form :

$$C=\left[\begin{array}{ccc}{C}_{11}& C_{12}& C_{16}\\ C_{12}& C_{22}& C_{26}\\ C_{16}& C_{26}& C_{66}\end{array}\right]$$

(3) Common symmetry classes

Material symmetry	Independent constants	Typical examples
Triclinic	21	General anisotropic crystal or composite
Monoclinic	13	Certain layered crystals and laminates
Orthotropic	9	Wood, woven composites, plywood
Transversely isotropic	5	Unidirectional fiber composites, graphite
Cubic	3	Silicon, copper, aluminum single crystals
Isotropic	2	Glass, isotropic polycrystals, random particulate composites
For the orthotropic material[2], An orthotropic material has three orthogonal symmetry planes. Since the symmetric plane yields $C=A_{\epsilon }^{T}C{A}_{\epsilon }$ (Derivation can be found in [2:1]). The stiffness at 4,5,6 out of diagonal is 0 :

$$\begin{array}{}\text{(1.3.1)}& {\mathbf{C}}_{\mathrm{ortho}}=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{12}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right].\end{array}$$

A transversely isotropic material^[3], with symmetry axis $3$, the matrix is :

$$\begin{array}{}\text{(1.3.2)}& {\mathbf{C}}_{\mathrm{TI}}=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{12}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{C_{11}-C_{12}}{2}}\end{array}\right]\end{array}$$

Cubic symmetric materials :

$$\begin{array}{}\text{(1.3.3)}& {\mathbf{C}}_{\mathrm{cubic}}=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{12}& 0& 0& 0\\ C_{12}& C_{11}& C_{12}& 0& 0& 0\\ C_{12}& C_{12}& C_{11}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{44}\end{array}\right]\end{array}$$

For the isotropic materials, we further have $C_{44}=\frac{C_{11}-C_{12}}{2}$, so we can express them to $\lambda$ and $\mu$ :

$$\begin{array}{}\text{(1.3.4)}& {\mathbf{C}}_{\mathrm{iso}}=\left[\begin{array}{cccccc}\lambda +2\mu & \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2\mu & 0& 0& 0\\ 0& 0& 0& \mu & 0& 0\\ 0& 0& 0& 0& \mu & 0\\ 0& 0& 0& 0& 0& \mu \end{array}\right]\end{array}$$

2. Transformation for Tensors and Voigt notations

(1) Engineering Strain vector

For strain with engineering strain notation, we should use :

$$\begin{array}{}\text{(2.1.1)}& {ϵ}_v=\left[\begin{array}{c}{ϵ}_{11}\\ {ϵ}_{22}\\ {ϵ}_{33}\\ 2{ϵ}_{12}\\ 2{ϵ}_{23}\\ 2{ϵ}_{13}\end{array}\right]\end{array}$$

This is because, the shear stress is ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{jk}$, when $j\ne k$, the component is :

$$\begin{array}{}\text{(2.1.2)}& C_{ijkl}{\epsilon }_{kl}+C_{ijlk}{\epsilon }_{lk}\end{array}$$

That is, when computing the stress, we should use the engineering strain vector to keep the expression consistent. So in the general Hook's Law^[1:1][4], we have :

$$\begin{array}{}\text{(2.1.3)}& \mathit{\sigma }=\mathbb{C}:{\mathit{\epsilon }}_v\end{array}$$

where ${\epsilon }_v$ is engineering strain vector.

(2) Implementation

We can transfer the format of them using the following two functions :

def tensor_to_voigt(tensor, *, engineering_shear: bool = False):
    r"""Return the Voigt vector of a symmetric tensor.

    The component order is ``[xx, yy, xy]`` in 2-D and
    ``[xx, yy, zz, xy, xz, yz]`` in 3-D. Set ``engineering_shear`` to
    ``True`` to multiply shear entries by two, as required for engineering
    strain vectors.
    """
    dim = tensor.ufl_shape[0]
    shear_scale = 2 if engineering_shear else 1
    if dim == 2:
        return ufl.as_vector(
            (tensor[0, 0], tensor[1, 1], shear_scale * tensor[0, 1])
        )
    if dim == 3:
        return ufl.as_vector(
            (
                tensor[0, 0],
                tensor[1, 1],
                tensor[2, 2],
                shear_scale * tensor[0, 1],
                shear_scale * tensor[0, 2],
                shear_scale * tensor[1, 2],
            )
        )
    raise ValueError("Voigt conversion supports only two- and three-dimensional tensors.")


def voigt_to_tensor(voigt, dim: int, *, engineering_shear: bool = False):
    r"""Return the symmetric tensor associated with a Voigt vector.

    The input component order is ``[xx, yy, xy]`` in 2-D and
    ``[xx, yy, zz, xy, xz, yz]`` in 3-D. Set ``engineering_shear`` to
    ``True`` when the input stores engineering shear entries, which are divided
    by two before constructing the tensor.
    """
    shear_scale = 0.5 if engineering_shear else 1
    if dim == 2:
        return ufl.as_matrix(
            ((voigt[0], shear_scale * voigt[2]),
             (shear_scale * voigt[2], voigt[1]))
        )
    if dim == 3:
        return ufl.as_matrix(
            (
                (voigt[0], shear_scale * voigt[3], shear_scale * voigt[4]),
                (shear_scale * voigt[3], voigt[1], shear_scale * voigt[5]),
                (shear_scale * voigt[4], shear_scale * voigt[5], voigt[2]),
            )
        )
    raise ValueError("Voigt conversion supports only two- and three-dimensional tensors.")

---

1. https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/06_LinearElasticity/06_Linear_Elasticity_03_Anisotropy.pdf ↩︎ ↩︎

2. https://en.wikipedia.org/wiki/Orthotropic_material ↩︎ ↩︎

3. https://en.wikipedia.org/wiki/Transverse_isotropy ↩︎

4. https://en.wikipedia.org/wiki/Hooke's_law ↩︎